function lu_experiment()
% LU_EXPERIMENT LU分解综合实验
% 
% 本实验涵盖LU分解的各种算法和应用

fprintf('=== LU分解综合实验 ===\n\n');

%% 实验1：基本LU分解
fprintf('实验1：基本LU分解\n');
fprintf('目标：理解LU分解的基本原理\n\n');

% 测试矩阵
A1 = [2, 1, 1; 4, 3, 3; 8, 7, 9];
fprintf('测试矩阵 A1 =\n');
disp(A1);

% 手工LU分解
[L1, U1] = lu_basic(A1);
fprintf('L矩阵 =\n'); disp(L1);
fprintf('U矩阵 =\n'); disp(U1);

% 验证分解
error1 = norm(L1 * U1 - A1, 'fro');
fprintf('分解误差 ||L*U - A|| = %e\n', error1);

% 与MATLAB内置函数比较
[L1_matlab, U1_matlab] = lu(A1);
fprintf('与MATLAB LU的差异: L差异 = %e, U差异 = %e\n', ...
        norm(L1 - L1_matlab, 'fro'), norm(U1 - U1_matlab, 'fro'));

%% 实验2：选主元的重要性
fprintf('\n实验2：选主元的重要性\n');
fprintf('目标：理解为什么需要选主元\n\n');

% 小主元矩阵
epsilon = 1e-10;
A2 = [epsilon, 1; 1, 1];
fprintf('小主元矩阵 (ε = %e):\n', epsilon);
disp(A2);

% 不选主元的LU分解
try
    [L2_no_pivot, U2_no_pivot] = lu_basic(A2);
    fprintf('不选主元成功\n');
    fprintf('L =\n'); disp(L2_no_pivot);
    fprintf('U =\n'); disp(U2_no_pivot);
    
    % 检查数值稳定性
    cond_L = cond(L2_no_pivot);
    cond_U = cond(U2_no_pivot);
    fprintf('L的条件数: %.2e, U的条件数: %.2e\n', cond_L, cond_U);
    
catch ME
    fprintf('不选主元失败: %s\n', ME.message);
end

% 部分选主元LU分解
[L2_pivot, U2_pivot, P2] = lu_partial_pivot(A2);
fprintf('部分选主元结果:\n');
fprintf('P =\n'); disp(P2);
fprintf('L =\n'); disp(L2_pivot);
fprintf('U =\n'); disp(U2_pivot);

% 验证 PA = LU
error2 = norm(P2 * A2 - L2_pivot * U2_pivot, 'fro');
fprintf('选主元分解误差: %e\n', error2);

%% 实验3：Cholesky分解
fprintf('\n实验3：Cholesky分解\n');
fprintf('目标：掌握对称正定矩阵的特殊分解\n\n');

% 构造对称正定矩阵
A3_temp = randn(4, 4);
A3 = A3_temp' * A3_temp + eye(4);  % 确保正定
fprintf('对称正定矩阵 A3 (%d×%d):\n', size(A3));
disp(A3);

% 检查正定性
eigenvals = eig(A3);
fprintf('特征值: [%.4f, %.4f, %.4f, %.4f]\n', eigenvals);
fprintf('最小特征值: %.4f (正定要求 > 0)\n', min(eigenvals));

% Cholesky分解
L3_chol = cholesky_decomp(A3);
fprintf('Cholesky因子 L =\n');
disp(L3_chol);

% 验证分解
error3 = norm(L3_chol * L3_chol' - A3, 'fro');
fprintf('Cholesky分解误差: %e\n', error3);

% 与LU分解比较运算量
n = size(A3, 1);
flops_cholesky = n^3/3;
flops_lu = 2*n^3/3;
fprintf('运算量比较:\n');
fprintf('Cholesky: %.0f, LU: %.0f, 节省: %.1f%%\n', ...
        flops_cholesky, flops_lu, (1-flops_cholesky/flops_lu)*100);

%% 实验4：带状矩阵LU分解
fprintf('\n实验4：带状矩阵LU分解\n');
fprintf('目标：利用矩阵结构提高效率\n\n');

% 三对角矩阵
n = 6;
A4 = diag(4*ones(n, 1)) + diag(-ones(n-1, 1), 1) + diag(-ones(n-1, 1), -1);
fprintf('三对角矩阵 A4 (%d×%d):\n', n, n);
disp(A4);

% Thomas算法（三对角LU分解）
a = -ones(n-1, 1);  % 下对角线
b = 4*ones(n, 1);   % 主对角线
c = -ones(n-1, 1);  % 上对角线

[L4_thomas, U4_thomas] = thomas_lu(a, b, c);
fprintf('Thomas算法结果:\n');
fprintf('L =\n'); disp(L4_thomas);
fprintf('U =\n'); disp(U4_thomas);

% 验证
error4 = norm(L4_thomas * U4_thomas - A4, 'fro');
fprintf('Thomas算法误差: %e\n', error4);

% 与通用LU分解比较
[L4_general, U4_general] = lu(A4);
fprintf('与通用LU分解的差异: %e\n', ...
        norm(L4_thomas * U4_thomas - L4_general * U4_general, 'fro'));

%% 实验5：数值稳定性分析
fprintf('\n实验5：数值稳定性分析\n');
fprintf('目标：分析不同方法的数值稳定性\n\n');

% 病态矩阵
A5 = hilb(5);  % Hilbert矩阵
fprintf('Hilbert矩阵 A5 (5×5):\n');
disp(A5);
fprintf('条件数: %.2e\n', cond(A5));

% 不同方法求解线性系统
b5 = ones(5, 1);
fprintf('求解 A*x = b，其中 b = ones(5,1)\n');

% 方法1：直接求解
x5_direct = A5 \ b5;

% 方法2：LU分解求解
[L5, U5, P5] = lu(A5);
x5_lu = U5 \ (L5 \ (P5 * b5));

% 方法3：Cholesky分解求解（Hilbert矩阵是正定的）
L5_chol = chol(A5, 'lower');
y5 = L5_chol \ b5;
x5_chol = L5_chol' \ y5;

% 比较解的精度
fprintf('解的比较:\n');
fprintf('直接求解: [%.6f, %.6f, %.6f, %.6f, %.6f]\n', x5_direct);
fprintf('LU分解:   [%.6f, %.6f, %.6f, %.6f, %.6f]\n', x5_lu);
fprintf('Cholesky: [%.6f, %.6f, %.6f, %.6f, %.6f]\n', x5_chol);

% 残差分析
residual_direct = norm(A5 * x5_direct - b5);
residual_lu = norm(A5 * x5_lu - b5);
residual_chol = norm(A5 * x5_chol - b5);

fprintf('残差比较:\n');
fprintf('直接求解: %e\n', residual_direct);
fprintf('LU分解:   %e\n', residual_lu);
fprintf('Cholesky: %e\n', residual_chol);

%% 实验6：性能测试
fprintf('\n实验6：性能测试\n');
fprintf('目标：比较不同方法的计算效率\n\n');

sizes = [100, 200, 500];
fprintf('矩阵规模    LU分解    Cholesky    Thomas\n');
fprintf('--------    ------    --------    ------\n');

for i = 1:length(sizes)
    n_test = sizes(i);
    
    % 生成测试矩阵
    A_dense = randn(n_test, n_test);
    A_dense = A_dense' * A_dense + eye(n_test);  % 对称正定
    
    % 三对角矩阵
    A_tri = diag(2*ones(n_test, 1)) + diag(-ones(n_test-1, 1), 1) + ...
            diag(-ones(n_test-1, 1), -1);
    
    % LU分解时间
    tic;
    [L_test, U_test] = lu(A_dense);
    time_lu = toc;
    
    % Cholesky分解时间
    tic;
    L_chol_test = chol(A_dense, 'lower');
    time_chol = toc;
    
    % Thomas算法时间
    a_test = -ones(n_test-1, 1);
    b_test = 2*ones(n_test, 1);
    c_test = -ones(n_test-1, 1);
    
    tic;
    [L_thomas_test, U_thomas_test] = thomas_lu(a_test, b_test, c_test);
    time_thomas = toc;
    
    fprintf('%6d      %6.4f    %8.4f    %6.4f\n', ...
            n_test, time_lu, time_chol, time_thomas);
end

%% 实验7：应用实例
fprintf('\n实验7：应用实例 - 电路分析\n');
fprintf('目标：使用LU分解分析电路网络\n\n');

% 简单电阻网络的节点分析
% 节点1: (V1-V2)/R1 + (V1-V3)/R2 = I1
% 节点2: (V2-V1)/R1 + (V2-V3)/R3 = 0
% 节点3: (V3-V1)/R2 + (V3-V2)/R3 = -I1

R1 = 1; R2 = 2; R3 = 3;  % 电阻值
I1 = 1;  % 电流源

% 构造节点导纳矩阵
G = [1/R1 + 1/R2, -1/R1, -1/R2;
     -1/R1, 1/R1 + 1/R3, -1/R3;
     -1/R2, -1/R3, 1/R2 + 1/R3];

I = [I1; 0; -I1];  % 电流向量

fprintf('节点导纳矩阵 G =\n');
disp(G);
fprintf('电流向量 I = [%g; %g; %g]\n', I);

% 使用LU分解求解节点电压
[L_circuit, U_circuit, P_circuit] = lu(G);
V = U_circuit \ (L_circuit \ (P_circuit * I));

fprintf('节点电压解: V = [%.4f; %.4f; %.4f] V\n', V);

% 计算支路电流
I_R1 = (V(1) - V(2)) / R1;
I_R2 = (V(1) - V(3)) / R2;
I_R3 = (V(2) - V(3)) / R3;

fprintf('支路电流:\n');
fprintf('I_R1 = %.4f A\n', I_R1);
fprintf('I_R2 = %.4f A\n', I_R2);
fprintf('I_R3 = %.4f A\n', I_R3);

% 验证基尔霍夫电流定律
KCL_check = [I_R1 + I_R2 - I1; -I_R1 + I_R3; -I_R2 - I_R3 + I1];
fprintf('KCL验证 (应该为零): [%.2e; %.2e; %.2e]\n', KCL_check);

%% 实验总结
fprintf('\n=== 实验总结 ===\n');
fprintf('1. LU分解是求解线性系统的基础方法\n');
fprintf('2. 选主元策略对数值稳定性至关重要\n');
fprintf('3. Cholesky分解适用于对称正定矩阵，效率更高\n');
fprintf('4. 利用矩阵结构（如带状）可显著提高效率\n');
fprintf('5. 不同方法在病态问题上表现不同\n');
fprintf('6. 实际应用中需要根据问题特点选择合适的方法\n\n');

end

% 辅助函数
function [L, U] = lu_basic(A)
% 基本LU分解（无选主元）
n = size(A, 1);
L = eye(n);
U = A;

for k = 1:n-1
    if abs(U(k, k)) < eps
        error('主元为零，需要选主元');
    end
    
    for i = k+1:n
        L(i, k) = U(i, k) / U(k, k);
        U(i, k:n) = U(i, k:n) - L(i, k) * U(k, k:n);
    end
end
end

function [L, U, P] = lu_partial_pivot(A)
% 部分选主元LU分解
n = size(A, 1);
L = eye(n);
U = A;
P = eye(n);

for k = 1:n-1
    [~, max_idx] = max(abs(U(k:n, k)));
    max_idx = max_idx + k - 1;
    
    if max_idx ~= k
        U([k, max_idx], :) = U([max_idx, k], :);
        P([k, max_idx], :) = P([max_idx, k], :);
        if k > 1
            L([k, max_idx], 1:k-1) = L([max_idx, k], 1:k-1);
        end
    end
    
    for i = k+1:n
        L(i, k) = U(i, k) / U(k, k);
        U(i, k:n) = U(i, k:n) - L(i, k) * U(k, k:n);
    end
end
end

function L = cholesky_decomp(A)
% Cholesky分解
n = size(A, 1);
L = zeros(n, n);

for k = 1:n
    sum_sq = sum(L(k, 1:k-1).^2);
    if A(k, k) - sum_sq <= 0
        error('矩阵不是正定的');
    end
    L(k, k) = sqrt(A(k, k) - sum_sq);
    
    for i = k+1:n
        sum_prod = sum(L(i, 1:k-1) .* L(k, 1:k-1));
        L(i, k) = (A(i, k) - sum_prod) / L(k, k);
    end
end
end

function [L, U] = thomas_lu(a, b, c)
% Thomas算法的LU分解形式
n = length(b);
L = eye(n);
U = zeros(n, n);

U(1, 1) = b(1);
if n > 1
    U(1, 2) = c(1);
end

for i = 2:n
    L(i, i-1) = a(i-1) / U(i-1, i-1);
    U(i, i) = b(i) - L(i, i-1) * U(i-1, i);
    if i < n
        U(i, i+1) = c(i);
    end
end
end